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In this note, following suggestions by Tao [2j, we extend the randomized algorithm 
for linear equations over prime fields by Raghavendra jlj to a randomized algorithm for 
linear equations over the reals. We also show that the algorithm can be parallelized to 
solve a system of linear equations Ax = b with a regular n x n matrix A in time 0(n 2 ), 
with probability one. Note that we do not assume that A is symmetric. 

Let m, n G N with m < n and consider an m x n matrix A e j^ mx ™ as we n as a 
right-hand side vector b e R m . There are many applications in which it is known in 
advance that A has full row rank, i.e. the system of linear equations Ax = b has at least 
one solution. We are interested in solving such a system in the sense that we want to 



*The author is indepted to Ian Hawke, School of Mathematics, University of Southampton, for pointing 
out an error in a previous version of this note. 
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construct a vector x G R n that fulfils these equations, given the knowledge that A has 
full row rank. This has, of course, important applications for the case m = n. For this 
problem, we consider the algorithm described below. In what follows, let A have the row 
vectors a 1; . . . , a m G R n , i.e. 




We consider a random vector £ : — > R n defined on some probability space 
(9, J 7 , P), for which the following holds. 

Assumption 1 For arbitrary a G R™, a^O, and (5 G R we have 

Prob (a T £ = /3) =0. 

In other words, the random vector is not biased towards particular affine subspaces of IR™. 
Examples for corresponding distributions include the case in which each coordinate (i = 
1, . . . ,n) is independently drawn from a Gaussian distribution on R, or from a uniform 
distribution over a certain interval, or in which £ is continuous uniformly distributed on 
the unit sphere {x G R" XlLi 3 ^ = !}■ F rom the assumption, it follows readily that 

Prob (f = x) = 0. 

for all x G R™, as Prob (f = x) < £^ =1 Prob ((e«) T £ = x;) = 0, where e« are the 
Cartesian unit vectors (i = l,...,n). For technical reasons, we will also assume that 
£(F) is a measurable set (in the usual sense of the natural Borel a- Algebra of R n ) for all 

We are now ready to state the main algorithm. 

1. Input: (A, b) with matrix A as row vectors ai, . . . , a m G R n , and a right-hand side 
b G R m . 

2. Let Vi, v 2 , . . . , v n+1 G R n denote identically independent distributed samples of the 
random variable £. 

3. for k — 1, . . . , m do 
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(a) Choose n+1 random pairs (h,ji), . . . , (i n +i, jn+i) with ig < jg for £ — 1, . . . , n+ 
1, and all pairs unequal to each other. 

(b) for £= 1, . . .,n+ 1 do 

i. x £ := rec(v if , v jv a. k ,b k ) 

(c) if one of the calls to rec stops with failure, then STOP with failure 

(d) Otherwise, set := for £ = 1, . . . , n 

4. Output: vi, . . . , v n+1 . 

This algorithm makes use of the subroutine rec ("recombination"), defined as follows: 

1. Input: (u, v, a, (3) with vectors u, v, a e IR™ and a real number 0. 

2. if a T (u - v) = then STOP with failure 

3. Otherwise, set 

/3-a T v 

t := — 

a T (u — v) 

and set z := tu + (1 — t)v. 

4. Output: z. 

In what follows, we will show the following. 

Theorem 1 Suppose A has full row rank and that Assumption^ holds. 

1. With probability one, the randomized algorithm described above stops after m steps 
with output ve e R™ ; I — 1, . . . , n + 1, such that Av^ = b holds for £ = 1, . . . , n + 1. 

2. With probability one, the run time of the algorithm is bounded by 0(n 2 m) floating 
point operations. 

From this, the following corollary immediately follows. 
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Corollary 1 Consider a regular matrix A G IR/ 1 *' 1 , a right-hand side b 6 R™ and suppose 
that Assumption U\ holds. Then, with probability one, the randomized algorithm above 
solves the linear system of equations Ax — b in 0(n 3 ) floating point operations. 

As it can be clearly seen, Step 3 of the algorithm can be fully parallelized. As each 
call to rec costs 0(n) flops, we arrive at the main result of this note. 

Corollary 2 Consider a regular matrix A G H nxn , a right-hand side h G IR™, and sup- 
pose that Assumption [I] holds. Then the randomized algorithm above can then be paral- 
lelized such that, with probability one, it solves the linear system of equations Ax = b in 
time 0{n 2 ). 

We start the analysis with a straightforward result. 

Lemma 1 Consider vectors u, v, a G R n with and a real number ft. Then, either 
a T (u — v) =0 or the subroutine rec returns a vector z = tu + (1 — t)v with a T z = (3. 

Proof: By construction. □ 

Next, we consider the first i iterations of the algorithm. 

Lemma 2 Let Assumption^ hold, let 1 < % < m and let a 1; . . . , aj be linearly indepen- 
dent. Then, the following holds. 

1. With probability one the algorithm has not stopped with failure in the first i iterations 
and the vectors vg = (£ — 1, . . . , n + 1), produced in step i of the algorithm, 
satisfy ajv^ = bj for j = 1, . . . , % and I — 1, . . . , n + 1. 

2. Let a G K,™ be an arbitrary vector with a^O and let j3 G R be arbitrary. Then, for 
all I = 1, . . . , n+ 1, a T v^ = (3 holds with probability zero, where = v® denote the 
iteration vectors of the algorithm after step i. 

Proof. We show both claims by induction. 

1. i = 1: claim 2 follows directly from Assumption [TJ Claim 1 follows from Lemma [H 
as A has full rank and a^Vj — Vj) = holds with probability zero for all i,j = 
l,...,n+l,i^j. 
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2. i — > i + 1 < m: We start with claim 1. Suppose that satisfy ajv^ = bj for 
j = 1, . . . , i and £ — 1, . . . , n + 1. Let (£i, £2) be a randomly chosen pair of indices 
with 1 < £1 < £ 2 < n + 1. Without loss of generality, assume £1 = 1 and £ 2 — 2. 
If rec(vi, v 2 , aj + i, returns a vector x without failure, then = b i+ \ by 

Lemma HJ Also, x is a convex combination of \\ and v 2 and thus fulfils ajx = bj 
for j = 1, . . . , i But, due to claim 2, the call rec(vi, v 2 , a$+i, &j+i) returns without 
failure with probability one. This shows claim 1. 

It remains to perform the inductive step for claim 2. As above, let us choose the 
pair of vectors v 1; v 2 without loss of generality. Due to the induction hypothesis, 
we have, with probability one, 

X = — rVi +11- —, T V 2 

^'(vi-vs) V ^'(vi-viOy 

= ^tTT -T (( b i - a7v 2 ) Vl - {h - a7 Vl )v 2 ) 

a i l V l ~ V 2j 

and therefore a T x = (3 if and only if 

[hi - a7v 2 )a T vi = (h - a.7vi)a T v 2 + /3aJ (vt - v 2 ) 

holds. Thus, 
Prob (a T x = 0) 

= Prob ((bi - a_7v 2 )a T vi = (b t - a 4 T vi)a T v 2 + /3aJ(vi - v 2 )) 

/oo 
Prob ((bi - a7v 2 )a T vi = ( and ( = (h - a7vi)a T v 2 + /3a7(vi - v 2 )) d( 
-00 

/oo 
Prob ((h - a7v 2 )a T Vl = C) + Prob (C = (h - slJv 1 ) S l t v 2 + /3a] "(v a - v 2 )) 
-00 

where the existence of the integrals are guaranteed as £ maps measurable sets on 
measurable sets, by assumption. But 

Prob (C = (h - a7 Vl )a T v 2 + /3a7(vi - v 2 )) 

/oo 
Prob (C = (&i - a7vi)a T v 2 + 0rj and 77 = a7(v x - v 2 )) dr] 
-oo 

/oo 
Prob (( = (h - a7vi)a T v 2 + f3ri) + Prob (77 = a7(vi - v 2 )) dr/ 
-oo 

/oo 
Prob (C - f3rj = (h - a7 Vl )a T v 2 ) dr] 
-oo 



< 
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for all C G IR, which shows 
Prob (a T x = 0) 

/oo /*oo 
Prob {{h - a7v 2 )a T Vl = C) + / Prob (C - Pv = (h - a7 Vl )a T v 2 ) drj d(. 
-oo J —oo 

Now, for C e R, C ^ 0, 

Prob ((h - a7v 2 )a T Vl = C) 

/oo 
Prob ((&< - ajv 2 )9 = ( and a T Vl = 0) d9 
■oo 

/oo 
Prob {{hi - a ? T v 2 )fl = C) + Prob (a T Vl = 9) d9 
-oo 

/oo 
Prob ((6j - ajv 2 )9 = C) d# 
-oo 

/0 ^oo 
Prob (6< - aT V2 = C/#) dfl + / Prob (6, - aTv 2 = C/0) 
oo JO 



d# 



and 

Prob — a7v 2 )a T vx = 0) = Prob (6^ = v 2 or a T V! = 0) 

< Prob (6j = a7v 2 ) + Prob (a T vi = 0) 
= 0. 

In a similar fashion, it can be shown that 

Prob ((-prj= {h - a7 Vl )a T v 2 ) = 

for all C, rj G R. As a consequence, a T x = (3 holds with probability zero. It is clear 
that the same analysis can be conducted for all other pairs of vectors v^, v^ 2 with 
h < i 2 . □ 



Lemma El invoked for % = m, shows part 1 of Theorem [U It remains to discuss the 
complexity of the algorithm. The f or-loop is over m steps, and each step involves three 
calls to rec. Executing rec costs two inner products of vectors in R n , two multiplications 
of vectors with scalars and one vector additon, i. e. the complexity of a call to rec is 0(n). 
These considerations show Part 2 of Theorem [TJ 

Some remarks are in order. 
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• It is clear that the algorithm also works for complex matrices A e C mxn and com- 
plex right-hand sides b e C n . Again, no symmetry assumption on A is necessary. 

• Some bookkeeping shows that the big-0 constant of the run time of the algorithm 
is ca. 15. While this appears large as compared to the big-0 constant of Gaussian 
elimination, 1/3, note that 15n 2 < n 3 /3 for n > 45. 

• The algorithm is optimal in the sense that its run time is of the same order as its 
input size (A, b). 

• The algorithm does not need to access the row vectors ai, . . . , a m directly; instead, 
it suffices to provide a routine that computes the action ajv of a row a, on an 
arbitrary vector vGR" 

• If the algorithm stops with failure in step k, then we have aju^ = ajv^ = ajw^ = bk 
for j — 1, . . . , k — 1, i. e. the algorithm provides at least solutions to a subset of the 
system of equations. 

• Stability issues: part of the stability of the algorithm rests on the size of quantities 
of the form l/(a^ +1 (u — v)). It is, at present, unclear how this quantity can be 
bounded away from zero. 

• In the exposition above, exactly n + 1 vectors are iteration vectors within the 
algorithm. We can, of course, use more than n + 1 vectors to iterate over, and 
choose in each step L > n + 1 pairs of vectors , v i2 from the current iterates 
to feed into rec. This increases the complexity of the algorithm from 0(n 2 m) 
to O(Lnm). However, choosing the right pairs of iterates v^v^ in an adaptive 
fashion, possibly discarding results whose norm is too large, might alleviate the 
stability issues mentioned above. 

• Another way that might be useful to stabilize the method at hand is to measure 
the degeneracy of a pair (vj, Vj) chosen in an iteration. If, say, ||vj — Vj|| is smaller 
than a certain threshold, the pair can either be discarded, or Vj can be replaced by 
Vj + c(vj — Uj) for a certain c > 1. A value of < c < 1 can be chosen if ||vj — Vj|| 
grows too large. 

• In Step 3 (a), it is not necessary to always choose n + 1 pairs of indices (and thus 
generate resp. update all of the vectors Vi, . . . , v n+ i. Indeed, after k steps of the 
main loop, all those vectors are in the /c-dimensional affine subspace defined by the 
first k equations ajx = bi (i = !,...,&) and will remain in this subspace for all 
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further iterations. Thus, after step k, only n + 1 — k pairs are needed to generate 
corresponding n + 1 — k new vectors. 
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